Skip to main content

Full text of "Minimal-Time Ship Routing"

See other formats


UNCLASSIFIED 



LIBRARY 

RESEARCH REPORTS 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA 93943-5002 



Technical Report 
distributed by 




Defense 
Technical 



Information 
Center 



am, 



/Acquiring Information- 
mparting Knowledge 



Cameron Station 
Alexandria, Virginia 22304-6145 



UNCLASSIFIED 



UNCLASSIFIED 



Policy on the Redistribution of DTIC-Supplied Information 

As a condition for obtaining DTIC services, all information received from DTIC 
that is not clearly marked for public release will be used only to bid or perform 
work under a U.S. Government contract or grant or for purposes specifically 
authorized by the U.S. Government agency that is sponsoring access. Further, 
the information will not be published for profit or in any manner offered for sale. 

Non-compliance may result in termination of access and a requirement to 
return all information obtained from DTIC. 

NOTICE 

We are pleased to supply this document in response to your request. 

The acquisition of technical reports, notes, memorandums, etc., is an active, ongoing 
program at the Defense Technical Information Center (DTIC) that depends, in part, 
on the efforts and interest of users and contributors. 

Therefore, if you know of the existence of any significant reports, etc., that are not in 
the DTIC collection, we would appreciate receiving copies or information related to 
their sources and availability. 

The appropriate regulations are Department of Defense Directive 3200.12, 
DoD Scientific and Technical Information Program; Department of Defense 
Directive 5230.24, Distribution Statements on Technical Documents; American 
National Standard Institute (ANSI) Standard Z39.18-1987, Scientific and Technical 
Reports- Organization,Preparation,and Production;Department of Defense 5200. 1-R, 
Information Security Program Regulation. 

Our Programs Management Branch, DTIC-OCP, will assist in resolving any 
questions you may have concerning documents to be submitted. Telephone numbers 
for that office are (703) 274-6847, or DSN 284-6847. The Reference Services 
Branch, DTIC-BCR, will assist in document identification, ordering and related 
questions. Telephone numbers for that office are (703) 274-7633 or DSN 284-7633. 

DO NOT RETURN THIS DOCUMENT TO DTIC 



EACH ACTIVITY IS RESPONSIBLE FOR DESTRUCTION OF THIS 
DOCUMENT ACCORDING TO APPLICABLE REGULATIONS. 



UNCLASSIFIED 



Lro/j6\ 



MINIMAL-TIME SHI? ROUTING 



by 

V. S. BLEICX and P. D. FAUULNE3 
Professor* of Math«?xat,i cs and Mechanics 



COPY ^ OF 



HARD COPY %.? f r 
MICROFICHE $. c ^O 



3Z f 

Research Paper No. 46 
UNITS) STATES NAVAL POSTGRADUATE SCHOOL 
Monterey, California 
August 1964 



Min&aal-Tiae Ship Routing 

V. E. BLIXCX end P. 0. PAUL**** 

U. 5. Revel Poetf r*4»*te School. Moaterey. Calif. 

ABSTRACT 

A rocoot theory of UBiMl-tiat ship routing through t lee-dependent 
oc«tn wave hoifht end diroctioo fields ie put to e numerical tool 
by using a series of eeaidaily analyse* furnished by tho U. S. Ravy 
Floot Numerical Veether Facility. Tho interpolations o>ad integrations 
required trt found to bo feasible. A roouao of tho thoory is given. 

1 . Introduction 
Heltiner, Heailton and 'Arnason (1962) gave a relaxation aothod 

solution to tho probloa of ainiael-tiee routing of ships through ocoa: 
viti hoight and direction fiolds dopondont on tho ship locatioa coor- 
dinate* only. Tho thoory has boon extended by Faulkaor (1963) to 
tho caso where tho v»»t hoight and direction depend on tiae also. 
Tho present paper confronts this tiae-dependent theory with actual 
vtYt height and direction analyses froa the files of the U. S. Navy 
Pleet Numerical Voather Facility, and reports on practical probleao 
which had to be solred in a test of the theory. 

2. Polar Telocity diagraa 

A basic ingredient of the theory is the polar diagraa of Fig. 1 , 
giving the ship Telocity t as a function of the angle I between the 
ship's heading and the wave direction. A diagraa of this kind aust 
be specified for each ware height H. The points L, N and N on the 
diagraa correspond to the ship speed v in head wares, v in beaa 
waves, and v ' in following waves. Eapirical curves for these three 
speeds as functions of wave height E are available la the pioneer 

- V- 



rmmrtrti -mnmrift iti iri in mi n ■■ ■■■ mh,m» ■■■..Mr mffciwni 








Pig. 1. Polar 



■w. ,.,_..., .,.—., ..^.^.^a. ^^-,. ■- ^^>^a»aL^a-^tt^AUa-M-aMkat 



: 




WAVE 



MOTION 



/ diagram. 



v~rk of Jascs i19V). Hit ?2sZ-:i2 »;.ip type c ..- -». i;.s«a in Fig. 2. 
ive been chosen for use here. They :.ave the appearance of being trcs 
of hvperboi^*, at itu*. approximately. T.. . * -*• confirmed wr.en * Utst 
•quart! analysis shov#d that all three 3f v , v ar.d v f can be re- 
presented closely as functions of H by the h> pt-rbol ic arc 

2 * 



C 1* C 2 H ~* ^ c 3* c 4 H ^* u :" c o' fc ~ c 3"^ 



(1) 



vhere c Q is the point cossoq to all t:.r*»e speeds *h«-n H-0. The other 
four constants are related to the asymptotes { c»*c . ) H*c . ♦ €-» and 
(c*-c , ) H*c . -c, . It vii decided to construct the polar velocity diagra 
by fitting an ellipse to the points L. .M and X. This resulted in semi 
principal axes a=(r ♦y-)/2, b^ar v /(v v ) , and a distance to the ec- 
centric pole given ty c*lv f -v )/2. Note t.\at tne pole is not to 
be construed as a focus of the ellipse. A further least squares aaaly 
sis shoved that the seaiaxes a and b are closely representabl e also 
by hyperbolic functions of the form of ( '. ) , but not c which oust be 
calculated as c = (y.-v, )/2. 
3. Coordinate system 

The ocean wave height and direction data of the semidaily Fleet 
Numerical Veather Facility analyses are presented in a south-polar 
stereographic projection of the nortnern hemisphere upon a plane pass 
ing through the circle of 60° North latitude. A rectangular coordinat 
system is set up in this projecting plane with the Ox and Oy axes 
parallel to the projections of the meridians of 10° and 100° East 
longitude respectively. A 62 by 62 grid is constructed using these 
axes with x=y=31 defining the projection of the North pole. The mesh 
distance between grid lines corresponds to a distance of 381 km at 
60° North latitude where the projection is true. The radius of the 
equator's projection is 31.205 mesh units. The map scale factor m, 

- 2 - 



BLANK PAGE 






> mttmetmMiriitfir i if'wmat^ammi^Autmm^ ji it iiifti i fimn i iitfMil* if iiiiiitaitoaM^MiMtiBttfcAfc. - 1 1 nritiwnil 



SHIP SPI 




Fig. 2. Ship speed in 1 




mt\*mt&**** f < ' i i'- "' ■••^** 



. —> . . 














■ 



;' 



;, beam and head vaves. 



TliftMmr "i ■■riWiaiii'itniiii iitMUftw'tiliiii i irai nij-inrn niii iftin ii'uan iimi iifi-jri'riiin i ir'i ; " - r'irl r Hi- ""^'^ — • ■""■ ■**•-• 



J 



defined as the ratio of a differential distance in the Oxy plane to 
the corresponding differential distance on tr.e earth's surface, is 

= = [973.75*U-31) 2 *(y-3O 2 j/l043.6 (2) 

Let x,y be the coordinates of a ship's projection in the Oxy plane 
at time t. Then the projected speed of the ship is 

V(x.y,t.e)=v(H.e)a (3) 

where v(H,9) is the actual geographical ship speed of the polar 
velocity diagram of Fig. 1, and where H(x,y,t) is obtained by inter- 
polation in the Fleet Numerical Weather Facility grid wave height 
data. Since the stereographic p*ojection is a conformal transforma- 
tion preserving angles and the*r senses, the angle S is the sace in 
the Oxy plane as on the earth's surface. 
4. Resume of the theory 

Fig. 3 shows a ship at the point (x,y,t) in the stereographic pro- 
jection plane on a route from fixed initial point A at t=0 to fixed 
teminal point B at t=T. Tne elliptical polar velocity diagram for 
V=nv is plotted at this point by interpolation in the . semidai ly wave 
height H(x,y,t) and wave direction K(x,y,t) grid values of the Fleet 
Numerical Veather Facility. The special type of interpolation re- 
quired is discussed in the Appendix. The direction of the ship's 
velocity vector V=ix"fjy is the control angle p, which is to be chosen 
at each point so as to minimize the transit time T from A to B. The 
equations of motion of the ship's projection in the Oxy plane are 

Qpjssx-V cosp=0, 5-=y-V smp=0 (4) 

where V*| V| sV(x,y , t ,p) . The problem of minimizing T is equivalent to 

the Lagrange cal ulus of variations problem of requiring the integral 

T 



I . J Q (WX^+u^dt (5) 



- 3 - 



■'■"" 



BLANK PAGE 



•**m<miii,m 1 1 ii«iMiiM«h— ». ■ nTaf,giiitmimitmmmdUitiUim*m*Mimti&*it*i^^ f i i 





Fig. 3 . Ship moti or 











1 



r mn i 



A 
i 

= 



Br t = T 




t 



4 
4 



1 



I 






eographic Oxy plane 



»-ni*iiMrtHfiiatte i mn ii n i t ■! 1 1 1 \mlamm imUmtJtUktt m i n , - , vn f r - : ' -^-^ -^-l 



,»<• "i raiiii M 



to be stationary, where X(t) and u(t) ir« contiou* of Lagrangian 
multipliers. Let the time at the fixed terminal point B be varied 
to T+&T. The vanishing first variation of 1 is 

6l-AT*[X6x«HA6y]J- f Q (g ) 5x«K? 4 6 r KP 5 6p)dt-C. (6' 

The coefficients of 6x,6y,6p in 61=0 give the Euler equations (7). 
(8) and (9) , consisting of the adjoint equations 

9^«X + (\ cosp-ni sinp)V «0, (7 

9 «u+(\ cosp+u sinp)V «0, (8 

and the scalar product control equation 

S> 5 =a.V =0, (9 

where the adjoint vector a«i\+ju, and where V sdV/dpei (Y cosp-Vsin 

«rj (V sinp+Vcosp) is the tangent vector to the polar velocity diagn 

of Pig. 3. Eq. (9) implies the orthogonality of and V as sbovn i; 

Fig. 3. Eq.(9) may be written also in the form 

p=arctanU/X)+arctan(V /V) . (10 

P 

The fixed end points A and B imply that 

dx(0)«6x(0)«0, dy(0)=<5y(0)^0 f (11 

dx(T)«(iAt+6x) T «0, dy(T)«(y,lt+6y) T =0. (12 

Use of £qs.(1l) and (12) makes the remaining terms of (6) proportio 

to wT , whose coefficient gives the scalar product transversal lty 

condition 

(^•Vj T =1>0 (13 

meaningful for sign only because of the homogeneity of (7) and (8) 

£q.(13) implies that the angle (q-p) between V and < is acute, as 

shown in Fig. 3. A further implication of (10) and (13) is that the 

quadrant of p is such that 

cosp-(XV-iiV )/AR f 3inp«(XV +uV)/ar, (14 



■•=■ " — 



vhere a.|a|.(x 2 >u 2 )^ and R»JV |»(V 2 *V 2 )*. 

Tht simultaneous numerical integration of (4), (7), (8) and (M) 
is carried out together vitb a Nevton-Rapheon iteration ai follovs: 
Let X^ ,u 1 and X-.WU b * tv0 li n#ariv independent solution* of the 
adjoint Eqs.(7) and (8) corresponding to the columns of the matrix 



E(t)-^1 X T| 



(15) 



vhere E(0)»I is the identity matrix. The X t u of (14) are taken as the 
linear combinations X=X. cosG+X -sinct and u=u. cosa-ni^sina. The vana- 
tion 6p is found by total differentiation of (10) to be 

6 pss R|E|6a/A 2 (V 2 +2V 2 -W ) (16) 

mm P PP 

vhere jEj is the determinant of E. Assume that a solution of the 

ship motion Eqs.(4) has been found, corresponding to (7) , (8) and 

(14) for some value of a, which falls short of the fixed end point 

B at t«T by the coordinate differences Ax(T) and Ay(T) . Using this 

solution and holding T fixed, find the variation of the vanishing 

matrix integral ~ 

J L¥ lt <P 2 ]S(t)dt«0. (17) 

Since the columns of E(t) satisfy the adjoint Eqs.(7) and (8), one 
obtains the 1x2 matrix equation 

[6x,6y] T E(T)=J [(V cosp-Vsinp) ,(V sinp+Vcosp) ]E6pdt. (18) 

Substitution from (14) and (16) into (18) gives 

[6x,6y] T =[-u,X] T J6a (19) 



where 



/"T 

3 - jmrjio *' |E| 2 dt/ A 3 (V 2 + 2V p 2 -W pp ) . (20) 



Nov vary the terminal time from T to T+AT and substitute 

[6x,6y] T= [Ax,Ay] T -[x,y] T AT (21) 

- 5 - 



into (19) to obtain the Nevton-Rapheon equations 

i(T)AT-Ju(T)6aȣx(T) 

(22) 
y(T)ATW\(T)6a«-ly(T) 

for the determination of AT and 6a on a varied trajectory vhich 
attempts to correct the errors Ax(T) and «ly(T) . The iteration to 
luccmiTt varied trajectories is continued until the terminal 
errors are acceptable. A suitable initial guess for the angle a is 
the inclination angle of the straight line from A to B. 
5. Nuaencal example 

Ten successive semidaily analyses of ware height and direction, 
starting at 062 on 4 May 1963, vere furnished by the Fleet Numerical 
Weather Facility. The maximum speed of the chosen P2-S2-R2 ship type 
is 19.6 knots. This combination of data precluded a trip of great 
length. It was decided to select an area of continued extreme wave 
height for the example. Such an area vas found centering at 30° North 
latitude and 162 East longitude. Figs. 4, 5 and 6 show tvo computed 
minimal-time ship tracks in the area, vith contours of vave height 
in feet and vave direction arrows. The arc traversed by the ship 
during the 6 hours preceding and/or the 6 hours following the time 
of each Figure is shown as a dashed curve. The minimal-time track 
AB required 2.488 days vith a 3.0^ saving over the geodesic track. 
The minimal-time track CD required 2.656 days vith a 1.3£ saving 
over the geodesic track. The severity of the sea conditions in the 
area preclude any more spectacular saving. The highly non-analytic 
nature of the vave height in the area vas found to affect the con- 
vergence of the Newton-Raphson iteration of (22). It vas found nec- 
essary to halve the values of AT and 6a in order to avoid a diver- 
gent oscillation. Resort to this stratagem vas found to be unneces- 

- 6 - 



BLANK PAGE 



.. ^..i .., ,-^J^ -■ , ...,.-■*■,--. ..-^- _.. 







Pig. 4. Sea conditions 

































• 



n 



».,^ lMIIM ii ■ \;,m ■ i -inn- ,ni - « *rtBmW0*+m 





nd 18Z on 4 May 1963. 



* i r > .»~m* ( t,+ sm , mdmimwrn • rif . M if fa ttt ia rrr n i1.ri . rto , - '■■HHftK l riiir ■ A W * lirt i i Hi i MMm ii T "T i n—ri d 




^-c-jL^^T 




Fig. 5. Sea conditions 



m.WM 



i jiUPJ Il l 'i il-.m. -.-liiU fi MHHP- ^ , 1 ^! ' Ujw. ' ^ ' "-A 1 H' l' 1 - 1 '- 1 ,** 1 ir " 



g i n ■■mfiri i ~— -■ * *.,.! n-fa., ■>>!■■* ■■!».- 



..<.■■*■,. 







id 18Z on 5 May 1963. 




■.■MWliVH- 



i r« . n i i « .« ,■ ». - ^ » m^ Mjhi^M^j 



. 



' 



* 




Pig. 6. Sea conditions 



■BH^ppw . i ii i .i i.i i-U P iw-.w."- ' u »- »j. " in. ^^w»^-^y;wM^pi^y»^p^^ _ : » »'i . ti.u»J .ui iiii i . . . 1 i«-),i"'. ' - wi 




5 



iiin**i wuriii 1 1 1 in,!.,.. ,1.1, ^ajwoMtatiiiatMfcrofcj.-, , «i"i».i , i ■■ 



',- 







P 



id 18Z on 6 May 1963. 



sary vhtn the vivi bright was nor* nearly analytic. 

6. Concluding remarks 

The numerical integrations involved in tne theory of minimal-time 
•hip routing through time-dtpendtnt ware fields arc found to be 
feasible. The necessary three-dimensional interpolations m the 
wave field data, discussed in the Appendix, present no problea. Con 
vergence problems may arise, but can be solved by the described de- 
layed approach to the limit. The authors can supply copies of their 
Fortran programs for ship routing and for the cubic-interpolation 
contouring of Figs. 4, 5 and 6. 

Acknowledgements . This work was supported by the Office of Naval 
Research. Aid from the U. S. Navy Fleet Numerical Veather Facility 
and from the Computer Facility, U. S. Naval Postgraduate School, 
is acknowledged. 

7. Appendix 

Some pertinent mathematical details are listed here. 
The geographic wave direction V, measured clockwise from the Nort 
must be converted to the unit vector l cosIC^j sinK in the stereo- 
graphic grid system by 

cosK«-[ (x-31 )cosv+(y-31 ) sinv]/r 
sinXa [ (x-31 ) sinv-(y-31 )cosv]/r 



(23) 



2 2 2 

where r »(x-31 ) +(y-31 ) . Then the derivative 

K =cosK(<JsinK/dx)-sinX(dcosK/dx) . 



(24) 



The derivatives V =mv +vm and V =av are obtained most conven- 

xxx P P 

iently by the implicit differentiation of the equation 

[v sin(p-K)/b] 2 +[(c+v cos( p-K) )/a] 2 =1 (25) 

of the elliptical polar velocity diagram, and noting that a,b,c 
are functions of H(x,y,t), and that K depends on x,y,t. The cod- 

- 7 - 



plexity of the result is reduced by introducing the parameter 6 

defined by 

sin9«b sin(q-K)/s* ▼ sin(p-K)/b 

(26) 

cos8-a cus(q-K)/s-[ v coi(p-K)-»c]A 

where s 2 *a 2 cos 2 (q-K)+b 2 sin 2 (q-K). 

The numerical integration of the adjoint Eqs.(7) and (8) demands 
an interpolation formula for H(x,y,t), cosK(x,y,t) and sinX(x,y,t) 
which guarantees the continuity of these functions and of their 
first space and time derivatives where any of x.y.t assume grid 
values. A 16-point interpolation formula to accomplish this is ob- 
tained from the 4x4 matrix F, whose four rows and columns of func- 
tion entries correspond to four successive x and y grid values re- 
spectively. The interpolation mesh cell is the central cell of the 
array, with x and y measured from the cell center, and with the mesh 
distance considered to be two units. The formula is 

F(x,y)=P(x)FPJ(y)/236 (27) 

where the matrix 

P(x)»[ ( 1-x) (x 2 -1 ) , (x-1 ) (3x 2 +2x-9) , (x+1 ) (9+2x-3x 2 ) , (xi-1 ) (x 2 -1 ) ] (28) 
and the prime indicates matrix transposition. Interpolation in the 
time dimension is accomplished by the similar 2-unit-mesh central- 
difference formula 

F(t)=[P(-3),F(-l),F(l),F())]?;(t)/l6 (29) 

which guarantees the continuity of F(t) and dF/dt at each end of the 
central time interpolation mesh. This formula is consistent with 
parabolic interpolation at the beginning or end of a time series, 
where central differences are not available. An interpolated vector 
i cosK(x,y ,t)t-j sin£(x,y,t) should be normalized before use. 

- 8 - 



REFERENCES 
Haltiner, G. J., Hamilton, H. D. and 'Arnason, G. # 1962: Minimal- 

tint ship routing. J. Applied Meteor ., 1, 1-7. 
Faulkner, F. D., 1963: Numerical methods for determining optimum 

ship routes. Navigation . 10, 351-367. 
James, B. V., 1957 (revised 1959): Application of vtvt forecasts 

to marine navigation . U. S. Nary Hydrographic Office. 



Fi 
Fi 
Fi 
Fi 
Fi 
Fi 



FIGURE LEGENDS 

1 • Polar Telocity diagram. 

2. Ship speed in following, beam and head waves. 

3. Ship notion in stsreographic Oxy plane. 

4. Sea conditions at 06Z and 1 8Z on 4 May 1963. 

5. Sea conditions at 06Z and 18Z on 5 May 1963. 

6. Sea conditions at 06Z and 1 8Z on 6 May 1963. 



- 9 - 



..JCB-BccICK - STAFF BOX 6 - 3 MINUTE EXPRESS - V7.5/6* 

PRCC*** ^INVCY 

EYVARS< U-Xl-l WARS < 2)*XMUl V VARS ( 3 ) *XLM2 YVARS(«»)*XM\j2 JJ«0 EW 

YVARS(S)*X YVAR<[6)-Y YVARS(7)*XJ YVARS<8)*S J>] *S 2 
01 "ENS I ON KHM9 v f2.10) 1 HT(19.12tlO).KK(19 v 12f 1C ) .CCSKI 1 9. 1 2 . 10). 3 

♦ $IMlM9.1?.10).fVARSm.OY(d). v C 18),C(*).TAUI306) t XI 300) .Y| 3001 9 % 

♦ CAPt A*( kO), PPM 00). CO I 3 00),*KliOO), VS(iOO) , »h ( 306 ) , X J ( J30 ) . * 

♦ S( 30C).AMfc v 6> 6 
LCUIVALESCE IKH.HT), (KK.COSK) 7 
«£A0 1, KH.KK.JJ, LHAX.KX$t,KYST.KXFN,KYFN.ALF ,T,FAC,FM\JL 8 

1 FOR-AT ( 10 ( 3912/ I • 60(381$/), 61 3* *F12.9 > 9 
PRINT 2. JJ, L*AX, KXST. KYST, KXFN, KYfK, MCt FMUL 

2 FCR»AT [|hO. 613, 2F12.9) 1 

nc u i- t . 19 1 

OELI • i-?a 1 

CO U J- t . 12 1 

DELJ ■ J*1 1 

ROOT » SwRTF(0ELl*0ELI ♦ CELJ*DELJ) 1 

CO U K« ■ , 1C 1 



MT { I t J.K1 • KM{ !• J V KI ] 

I.J.K) • 1 
ANCLE ■ ANCLE/57,2957 7951 



ANGLE • KK( I.J.K) • 10 2 



CCS - COSF(ANCLE) 

SIN ■ SINF(ANGLE) 

COSK{ I.J.K) • -OELI'COS ♦ OELJ»SIN)/RCOT 

SINKII,J,K) • C0ELI*SIN - OELJ*COS)/ROOT 

C( 1 ) ■ 0.0 

C(2) - 0.5 

C(3) -0.5 

C(U) • 1.0 2 

X( 1 ) • KXST 2 

Y( 1 ) ■ KYST 

TAU( 1 )« 0.0 

XFIN * KXFN 

YFIN • KYFN 

WH( 1 ) » HT(KXST,KYST. 1) 

CALL POLAR (COSK( KXST.KYST.l) , SINK ( KXST.KYST, 1) , WK ( 1 ) I 

XSTEP « T/FAC 

CALL ABC (WH( 1),AI,BI,CI,0AI,DBI,DCI) 

CAPLA^( 1 ) » 1.0 * 

S( 1 ) ■ 0.0 »» 

XJ( 1)* 0.0 * 

00 18 L-l.LMAX h 

- ALF k 

Q0( 1) - C • 57.2957 7951 * 

COSA « CCSF(ALF) % 

SINA m SINF(ALF) U 

COSQ ■ CCSA h 

SINQ * SINA * * 

COSQMK- C0SQ»C0SK(KXST,KYST,1 ) ♦ SI NQ»Sl NK < KXST, KYST, 1 ) U 

SINQ.MK* SIN3»C0S<(KXST,KYST,1 ) - COSQ*SINK( KXST, KYST, 1 ) 

A8S ■ AI«COSQMK 

ORO ■ 9I«SINQMK 5 

HYP » SCRTFt ABS«ABS ♦ OR0*OR0) 5 

SINB* ORC/HYP 

COSB* ABS/HYP 

VMAJ- AI«COSB - CI 

V M I N» P I • ^ I H fl 

VSt 1)«SQRTF(VMAJ»VMAJ ♦ VMIN«VMIN) 

CALL POLAR (VMAJ, VMIN.PMK) 

PP( 1 )• PMK ♦ WK(1 ) 6 

XVAR - 0.0 6 

00 5 1-1,8 6 

YVARSd ) » 0.0 6 

YVARS(l) • 1.0 t 

YVARS(U) « 1.0 t 

YVARSI5) « X( 1) t 

YVARS(6) ■ Y( 1) 6 

Nl - T/XSTEP ♦ 1.0 4 

XN1 * Nl t 

STEP - T/XN1 7 

N2 ■ Nl ♦ 1 7 

00 1U K«2,N2 7 

DO 7 1*1.4 7 

XC * XVAR ♦ C(I )*STEP 7 

DO 6 J«l,8 7 

YC(J)- YVARS(J) ♦ C» I )«AK(I-1,J) 7 

XLAM - YCID^COSA ♦ YC(3)»SINA 7 

XMU * YC(2)*C0SA ♦ YC(U)»SINA 7 

CLAM • SCRTF( XLAM»XLAM ♦ XMU«XMU) 7 

CALL TERP <HT,YC< 5 ) , YCI 6 ) , XC, M, HX.HY) i 



8 



80 
81 



CALL 
CALL 

CA' ' 

ROOT 

CK 

SK 

OKX 

OKY 

cosq 

INQ 
OSQM 
SINQH 
ABS 
CRD 
HYP 
SIN8 
COSB 
VMAJ 
VHIN 
V 
COSP 
SINP 
COST 
VBR 
OELX 
DELY 
EMFI 
CAPV 
OY(S) 
0Y16) 
EMFIX 
EMFIY 
0Y(8) 
B2HA2 
AMCCB 
RAT 
VP 
CAPVP 
OIV-1 

OCT 
RADIX 
0Y(7) 
FNUM 
DBDH 
RATX 
RATY 
VBRX 
VBRY 
CAPVX 
CAPVY 
0Y< 1) 
0Y(2> 
0Y( 3) 
DY(U) 
00 7 
AK( I , 
DO 8 



APC (h.pa.FB.F 
TERP itll* 9 Vti 
T£«»o , SINK,YC(5).YC(6), 

s;rtfich»;k ♦ sk»sm 

CK/RC3T 



v 0A.n8.0C) 

) .YC(6)#XC,C«t.C>CX t CHY) 
XC,SK,SK*,5KY) 



SK/RCOT 

C*»SKX - SK«CKX 



♦ SINQ«SK 

- CCSG«SK 



♦ ORO*ORO) 



CK«SKY - S*«C*Y 
XLA-/CLAH 
X^U/CL I* 

• COSC«CK 

• SINQ«CK 
FA«CCSQPK 
FC«SINCMK 
SwRTFI ABS* AeS 
OaC/HVP 
ADS/HYP 
FA.CCSB - FC 
FP«SINR 

SQRTF( VMAJ»VNAJ ♦ VMIN«VMIN) 
(CK*VMAJ - SK»V*IN)/V 
(S*»VMAJ ♦ CK»VMIN)/V 
V^AJ/V 

V/1R.7C2 181 818 
YC(S) - 26.0 
YC(6) ♦ 1.0 

(973. 75 ♦ OELX»DELX ♦ OELY^OEL Y) /l 0M3.638 7U3 
veR»E*FI 
CAPV^COSP 
CAPV«SINP 
0CLX/S21.819 
0ELY/S21.319 
CAPV 

cp«FB - FA»FA 
FA - FC»C3SB 
SINB»(FA»FC ♦ 
RAT • V 
RAT • CAPV 

♦RAT«»2-( V/AMCCB)»»2»( (VP»CCST/SINB)-(82MA2»SIhB*»2/FB))/FB 
YC( 1 )»YC(U) - YC(2)«YC(3) 
SCRTFl 1.0 ♦ RAT»RAT) 

(CAPV»DET«0ET/0IV) • ( R AOI X/CL AM ) ««3 
-RAT»APCC8 

♦ FA«08»SINB««2/FB - DC^COSB 

♦ HX*OBOH)/AMCCB 

♦ HY«OBDH)/AMCCB 



3715 
3715 



B2HA2«C0SB)/(F8«AMCCB) 



VBR»EMFIX 

VBR'EMFIY 
1MC0SP ♦ YC(21«SINP) 
YC(2)«SINP) 
YC(U)»SINP) 
YC(»»)«SINP) 



1)*C0SP 
3)»C0SP 
3)«C0SP 



♦ 



0A»C0S8»»2 
(FNUM»OKX 
(FNUN«UKY ♦ 
RATX»VBR 
RATY»VBR 
V9RX«E.rFI ♦ 
VBRY«E*FI ♦ 
-CAPVX»(YC ( 
-CAPVY»( YC( 
-CAPVX#( YC ( 

-capvymy: ( 

J-l ,8 
J) » STEP • DY(J) 
J-l, 3 

YVARS(J) -YVARS(J) ♦ { 

XVAR « XVAR ♦ STEP 

TAU(K) • TAU(K-I) ♦ STEP 

X(K) » YVARS( 5) 

YIK) - YVARS16) 

IF IN2-K) 80.81,80 

IF (LMAX-L) 82,81,82 

XLAH ■ YVARS( 1 NCOSA ♦ 
XMU ■ YVARS(2)»C0SA ♦ 

CLAM ■ SORTF (XLAM«XLAM 

CAPLAM(K) * CLAM 

CALL TERP (HT,X(X) f Y{K).TAU(K) ,WH(K),HX,HY) 
ABC (WH(K),FA, FB v FC.0A.0B,0C) 
TERP (COSK,X(K),Y(K),TAU(K) t CK,CKX,CKY) 
TERP ( SINK,X (K),Y(K) , TAUIK), SK,SKX,SKY) 
• SCRTF(CK<K ♦ SK«SK) 
CK/RCOT 



AKtl t J)*2.»AK{2,J)*2.»AH(3,Jl*AK(U,J) )/6. 



YVARS(3)«SINA 
YVARS(U)»SINA 
♦ XMU«XMU) 



CALL 
CALL 
CALL 

ROOT - 
CK • 
SK • 

COSQ « 

SINQ • 

COSQMK 

SINQMK 

ABS - 

ORO ■ 



SK/RCOT 
XLAM/CLAM 
XMU/CLAM 

- COSO»CK 

- SINOCK 
FA»COSCMK 
FB«SINQMK 



♦ SINQ»SK 
- COSQ«SK 



- 11 - 



92 
83 
8U 
85 
66 
87 
88 
89 
90 
91 
92 
11 
12 

13 

lu 
15 
16 

17 



HYP 
SIN8 

cose 

VMAJ 
VMIN 
VS<* 
CCSP 
SINP 
CALL 
CALL 
CALL 
XJ(K 

IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
N2 
T 



• SC 

• 0* 

• AP 

• FA 

• FP 
■ Si 

• (C 
» (S 
PCLA 
PCLA 

prLA 

• Y 

• Y 
(N2-K) 
( JJ) 
(XFIN- 
(X(K)- 
(X (K)- 
(Y(K)- 
( 1 CO- 
(YFIN- 
(Y(K)- 
(Y(K)- 
(X(K)- 
( 16.0- 

■ K 

■ TAU( 
PRINT 13. 
FORMM (1 
GO TO 81 
CONTINUE 
PRINT 16 
FORMAT { 1 
PRINT 
«ORMA 
D£LX 



A8S ♦ ORO«ORO) 



FC 



V*!N»VMIN) 
)/VS(K) 

)/VSIK) 

K)) 

KM 



•V*AJ ♦ 

CK»VMIN 
SINP,PP( 

SINQ,QQ( 

!yVARS( 1 )«YVARSIU) - YVARS(2)»YVARS( 3) ) 

82 



I5,21,8w 
16, 11,11 
) 11,11,86 
) 11,11,87 
) 1 1,11, U 
) 90,21,89 
) 91,11,11 
) 1 1.11,91 
) 11,11,92 
) 1 1,11, U 

EARLY N2» 12/1 



17 



X( 
DELY • Y( 
EMFI * (9 
CAPV » VS 
XOOT « CA 
YDOT ■ CA 
FINJ « XJ 
DETER- FI 
OIFX ■ XF 
OIFY • YF 
DIFT » FI 
OIFA « (X 
T « T 
ALF » AL 
PRINT 16 
PRINT 17, 

18 CONTINUE 
PRINT 19 

19 FORMAT( 1» 
♦8X2HPP8X2 

PRINT 20, 
♦ 

20 FORMAT (1 

21 STOP 
ENO 
SUPROLTIN 



H06X 
N2, 
10. 
N2) 
N?) 
73.7 
(N2) 
PV«C 
?WS 
(N2) 

NJ«( 
IN - 
IN - 
NJM 
DOT* 
♦ FM 
f ♦ 



2HN2 
ALF. 
4F15 
- 26 

♦ 1 
5 ♦ 

• EMF 

OSP 
tNP 

XDCT 
X(N 
Y(N 

XLAM 
UIFY 

ULO 
FMUL 



12X3HALF 

T.X(N2), 

.9 ) 

.0 

.0 

OELX»CEL 

1/18.702 



1UX1HT1UX1HX13X1HY) 
Y<N2) 



X ♦ DELY»0ELY)/10U3 
181 318 



638 7*3 



•XLAH ♦ YOOT«XMU) 

2) 

2) 

•OIFX ♦ 

- YOCT* 
IFT 
• OIFA 



XMU»0IFY)/CE1ER 
0IFX1/DETER 



N2,ALF, T,X(N2),Y(N2) 



CUX3HTAJ8X1HX9X1HY6X6HCAPLAM5X2HWH8X2HWK8X2HVS9X1HS8X2H 

HCQ/ ) 

(TAU( I ). X( I) ,Y( I),CAPLAM( I),WH( I),WK( I ), 

VS( I ).S( I ),XJ( I ),PP( I ),QQ( I ), I»1 ,N2) 
1F10.4) 



RA01 

A 

RA02 

8 

RA03 

RA04 

C 

OA 

08 

OC 



SC 

-c 

SC 
-0 
SC 
SC 
C 
-0 
-0 

c 



E ABC M 
RTF( (0.0 



120 
RTF( 
.126 
RTF( 
RTF( 
.125 
.120 
.126 
.125 



197 

(0.0 

879 

(0.0 

(0.0 

U2»» 

197 

879 

U2U 



, A, 8. C,0 

38 635 
666»H 

32 937 
519«H 

10 Ui» 1 

13 051 
U93«H 
666 -( 
519 -( 



ENO 
SUBROUTIN 



IF (X) 
IF (Y) 

P- O.C 
RETURN 
P«90.C 
RETURN 
P« -9C.0 



E POLAR 
0,3, 10 
8,U,6 



*93 -(0 
(X,Y,P) 



A, 08, DC) 
6U66»H-0 

20.7U8 
638U»H-0 

20.633 
95i»M-0. 
9lO»H-0. 

0.29} 
.038 635 
.032 937 

:§1S ii\ 



.492 U35 972 

910 236 - RA 

.U39 806 791 

U62 763 - RA 

066 47<? 3891 

193 791 63U 

**U4 e925 - R 

06«*66«ri - 

3638U»H - 

8956Q«M - C 

79100»H - 



l\ 



•H+2.099 340 872) 



122 9UC) 



9)«H*1.778 

D2 

)»m>0.371 204 63** 1 ) 

)«H«-0.814 886 3262) 

AD3 ♦ RAOU 

.2«»6 217 98635J/RAD 

.219 903 395951/RAC 



:8« 11? 



17V 

173 
173, 
174 
174 
174. 
174. 
17w.i 
175 
176 
177 
178 
179 
180 
181 
182 
183 
184 
185 
166 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
2C0 
201 
202 
XJ203 
204 
205 
206 
207 
203 
209 

1 
2 
3 

4 

5 

6 
7 
8 
9 



- 12 - 



is 

1 

I 

4 

5 
6 

7 
8 



9 

IJ 

5 

16 



RETUR 
P- 57 
IF (X 
IF (V 
p. p 

RETUR 
P- P 
END 

SU8R0 
OINEN 



TO 



(L 

id 



TO 



11 

L 

IF 

IF 

LL 

L 

CO 

LL 

L • 

M • 

N ■ 

XX • 

VY ■ 

TT • 

XP1 ■ 

XM1« 

YP1» 

YM1» 

TP1- 

TM« 

X2H« 

y2h« 
T2M« 
P( 1) 
P(21 
P(31 

P(U ) 
01 1 ) 
0(2) 
Q(3) 
Q(u) 
PX(U) 
PX( 1) 
PX(2) 
PX( 3) 

OY(U) 
QY{ 1) 
QY(2) 
0Y(3» 
DC 6 
00 5 
DO 5 

DO 16 
S( I 1 
D( I ) 
OROJK 
OROX( 
OROY( 
IF (L 
G - 

H « 

OUT » 
OUTX» 
OUTY* 
RETUR 
TM3 ■ 
OUT - 
OUTX» 
OUTY* 
RETUR 
TP3 • 
OUT « 
CUTX- 
OUTY- 
ENO 

EN 
> 10 1 



.2957 7951 • 
) I2.U.U 
) 15,13,13 

160. 

160.0 

TINE TERP ( 

ION FUNCM9 

ORDAIM 

INTFCT1 
ItO 
L) 3 2,3 

-2 

k 



ATANFC Y/X) 



FLiNC, X, Y, T, OUT. CUTX, CUTYl 
,17JO) T PTU.U).P(%) f C(»»I,PX(H),QY{i») ,OROU)t 

,or6y(*),s<*).oU) 



-i 

INTPt 
INTF( 
.0»(X 
.OMY 
.0»(T 
X ♦ 1 



P1»X* 

P1»VM 

P1»TM 

-XKl 

1(3. 

-2.» 

XP1 
-YMl 

((3. 

-2.* 

YP1 

(3.« 

u.» 

(9.« 
-U.« 
11. • 

(<=.• 
-u.» 

K«l, 

!•!• 

J-1, 
J) - 

1-1. 

• P( 1 

C(1 

«(D( 

•(D( 

■(Si 

L-21 

(3.0» 

G - U 

(TM1«( 

(TM1«( 

(TM1-I 

N 

TT - 

( TM3 

( t*3 

( TM3 
N 
TT ♦ 
( TP3 
( TP 
( TP 

C 
C 10 



X) - 2 

Y) - 2 

-INVF( 

-INTF( 

-INTF( 

.0 

.0 

.0 

.0 

.0 

.0 

1 
1 
1 

• X2M 
•xx*2. 
xx«xx 

• X2H 

• Y2* 

• YY*2. 
YY«YY 

• Y2H 
XX-1.) 
XX - P 
XX-11. 
XX - P 
YY-1.1 
YY - Q 
YY-11. 
YY - 



X)) - 1.0 
Y11 - 1.0 
T)l - 1.0 



) 

K) 

Kl 



UNC(M-f 
H 

)«PT( i 
UPTU 
1 )• P( 
11»PX( 
1 )«0Y( 

8 t 7 l 9 
TTO.O 

.0»TT 

G»ORD( 

G^ORDX 

G«ORDY 

3.0 
•( TM1» 
•{TM1» 
•(TM1» 

3.0 
•(TP1» 
•{ TP1» 
• ITPH 



l»XX-9. )*XH1 

♦ 18. - P(2) 

)»YY-9.1«Y*1 

♦ 18. - Q(21 

• XP1 
X(U) 
)»XP1 
X(2) 

• YP1 
YU) 
)»YP1 
Y(2) 



I,N*J,L-*M 



,h* 

1 ) ♦ 
1) ♦ 
1) ♦ 



P(2)»PT(2,I )♦ P 



0(2)*PT( I, 
D(2)» P(2l 
D<2)«PX(2> 
S(2)«QY(2) 




3)»PTC3tI)* 
3)»PT( I, 3h 
)• P(3) 
)»PX(3) 
1»QY(3) 



P(«)»PT(U, I) 
C(U)«PT( I ,u ) 

♦ 0( »»)• P(U) 1/256 

♦ D(U)«PX(U) )/12e 

♦ Sm«CYU) )/ 128 



1»TT - 9.0 



2)-T2*»ORD( 1 )1-TP1»(H«0R0( 3)-T2*^CR0f Wlll/16. 
(2>-T2*«0R0X( 1 11-TP1 •(H»CRDX(3)-T2M«CR0X(U) 11/16. 
12)-T2H»ORDY(1))-TP1«(H«GROY(3>-T2M»CROY(»») 11/ 16. 



oroi ii-; 
oroxc u-: 
oroy(I)-; 



.0«TPUORD( 211 
!.0*TP1«0R0X<2> 1 
!.0«TPl«OROY(2l 1 



0R0( M-2.0»TMW0R0l 

0R0X(*)-2.0«TM1«0R0X( 

ORDYU)-2.0»TM1«OROY( 



♦ T2*»0R0( 311/8.0 

♦ T2W«GRDX( 3) 1/8.0 

♦ T2*«CR0Y{3!)/8.0 



♦ T2^»CRD( 21 1/8.0 

♦ T2*«CR0X(2> 1/8. C 

♦ T2*«ORDY(2) 1/8.0 



5-1. 362675810 %. 9770*031 150. OOOOOOCOO 0.500000CCO 



PROGRAM STRAIGHT 



- 13 - 



C TO CCNVPBT PQTGaA* MNVCY TO PROGRAM STRAIGHT 

C (1) REPLACE wisvnv CARDS (U3-55) BY FOLLOWING CAROS 

C COST --Sl?.K(KXST f KVST, 1 ) % 

C OR COST ■ CCS«(KXST, KYST, 1) % 

C SINT — CCSMKxST,KYST,l ) % 

C OR SINT — SINK (KXST.KYST. 1) fc 

C CUAO • 8I«RI«C0ST»CCST ♦ AI • A I «S I NT«S INT % 

C RAOI • SwRTF(CUAD - ( C I »S INT ) «»2 ) * 

C SINB ■ SINT«( AJ.RA0I-BI»CI»C0ST)/GUA0 «» 

I COSP - (flt»COSTM ACI*Al*Cl«SINT«SlNT)/CUAO *» 

C 12) REPLACE *INVOY CAROS (89-97) BY FOLLOWING CARCS 

C OR COST I Ik 8 

C SINT «-CK 9 

C OR SINT — S* 9 

C CUAO • FB«FR»COST«COST ♦ FA»FA« SINT«S INT 9 

C RAOI ■ SwRTFI 0UA3-(FC»SINT)«»2) 9 

C SINB ■ SINT«(FA«RACI - F 8 «FC 'COST ) /QU AO 9 

C COSP « (FB»COST«RADI ♦ F A»FC» SINT«S INT ) /GUAO 9 

C (3) REPLACE NiNVOY CAROS (153-161) BY FOLLOWING CAROS 

C COST «-SK 

C OR COST ■ CK 

C SINT «-CK 

C OR SINT »-SK 

C QUAO * FP»FB»CCST#CCST ♦ FA»F A»SI NT»S INT 

C RAOI • SCRTF( CUAO-(FC»SINT)««2) 

C SINB - SINT#( FA»RAOI - F B»FC»COST ) /QU AO 

C COS8 « (f B»COST»RACI ♦ F A»FC»S INT»S INT ) /QLiAO 

C (U) REPLACE MINVOY CAROS (197-198) BY FOLLOWING CAROS 



C FF - (5.CO-Y(N2-l ) )/( Y(N2)-Y(N2-1 )) 
C OR FF • ( 13.0-X(N2-1 ) )/(X(N2)-X(N2-l ) ) 
C T - FF»TAU(N2) ♦ ( 1.0-FF)»TAU(N2-1) 



PROGRAM CONTOUR Q 

DIMENSION ABS(900 ) .ORO (9CO) , KH( l9,12),HT(19,12),KK(19,12),PT(12) t 
♦COSK( 19,12)tSINK(19.12),RX(U8),RY(«»8),C0N(U f u) t IT(12),THl2),C(3) 2 
EQUIVALENCE (IT.TI), I L A , AL ) 3 

READ 1, GS, IHMAX, IDELH. I XMIN , I XMAX. I YMI N, I YNAX ,KH, KK, T I , AL * 

1 FORMAT i F±.2i 613/, 12(3812/), 2(l0A8/> > 5 
00 2 J- 1 ,12 6 
0021-1,19 7i 

2 MT( I. J) • KHI I, J) 8 
PRINT 3, OS. IHMAX, IDELH, IXMIN, IXMAX, IYMIN, IYMAX 9 

3 FORMAT ( IhO. FU.2. 613 ) 
READ 26, (ABS(I), ORO(I), 1*1,5) 

26 FORMAT ( 10F3. 1) 

CALL CRAM (5, ABS, ORO, 1 ,0, LA, I T, 1 . , 1 . , 0,0, 2,2,8, 8,0, LAST) 

C » DS»DS 

DO 22 IY- IYMIN, IYMAX 

Y - IY 
00 22 IX- IXMIN, IXMAX 

CALL COEF (HT, IX, IY. PT, CON) 

DO 22 IH- 3, IHMAX, lOELH 
H * IH 

NR — C 

CALL BOOT (CON, 1, H, 0, KER) 

IF (KER) 8.8,% 
H 00 7 J-1,KER 

ABS( 1 ) - X ♦ (0(J)*1. 01/2.0 

ORD( 1 ) » Y 

IF (NR) 33,35,33 
33 DO 3U I-l.NR 

A - AeS( 1 ) - RXU ) 

B - CRC( 1 ) - RYU ) 

IF (A«A ♦ B«B - C) 7,7, 3H 
3U CONTINUE 

35 CALL GRAC ( PT , ABS ( 1 ) .ORO ( 1 ) ,0X, DY,OUAO,OELX .OELY.OS ,H,0) 
IF (QUAO - l.E-10) 7,7,5 

5 IF (DELY) 6,36,36 

6 OS - -OS 
OELX- CELX - DX - OX 
DELY- DELY - DY - DY 
IF (OELY) 7,36,36 

36 ABS(2) - ABS(T) ♦ OELX 
0RD(2 ) - ORQ( 1 ) ♦ OELY 
CALL CUT (NK, ABS,ORO,RX,RY,OS,PT,X,V, H) 

7 CONTINUE 
IF (KER) 13.13,9 

8 CALL ROOT (CON, 2, H, 0, KER) 

9 00 12 J-1,KER 



37 



ABSC 1 1 • X 5 

CR0(1) • Y ♦ ( 0( J ) ♦ 1 . 01/2.0 5 

IF (NR) 37,39,37 5 

00 38 I-l.NR 5 

A • ABS( 1 ) - RX(l ) 5 



p • 0«0( 1 ) - RYU ) 

|F (A-A ♦ B*B - t> 12,12,38 

8 CONTINUE 

9 CALL CRAO (PT.ABSt 11.0*0(1). OX, DY, QUAO, OELX, OELY, OS, H,0) 
IF (OLAC - ) t E-0) i2,i2,fo 
IF (OELX) 11,U0,%0 
OS • -OS 

8 Ely* DFLY - dy - oy 
elx- Ctix - ox - ox 
IF (OELX) 12.uo.uc 

*0 ABS(2 ) - A8S( n ♦ OELX 
0R0(2 ) • 0R0( 1) ♦ OELY 



!? 



CALL CUT (NR,ABS,OR0,RX,RY,OS,PT,X,Y,H) 
. 12 CONTINUE 

13 CALL ROOT (CON. 3, M, 0, KEH) 

IF (KER) 18,18,1% 
1% 00 17 J-1,K£R 

APSC1 ) • X ♦ (0(J )«1. 01/2.0 

ORD( 1 > • Y ♦ 1.0 

IF (NR) Wl.U3.Hl 
fcl 00 *2 I-l.NR 

A • A6S( ) ) - RX(I 1 

8 • OROI 11 - RYU ) 

IF (A«A ♦ 6*8 - CI 17,17,%2 



*2 CONTINUE 

U3 CALL CRAC (PT.ABSt 1),CRQ( U.OX.OY, QUAO, OELX, OELY, OS, h,0) 
IF (QUAO - 1. £-10 1 17,17,1$ 

16 



l.E-10) 
IF (OELY) UU.UU, 16 « 

OS • -OS « 

OELX- OELX - OX - OX « 

OELY- CELY - OY - DY « 

IF (OELY) UU.UU, 17 ■ 

M ABS«2) ■ ABS( 1) ♦ OELX 
0R0(2 ) • CROC 1) ♦ OELY 

CALL CUT (NR, A8S,OR0,AX,RY,0S,PT,X,Y,M) 
CONTINUE 

KER) 



Is 7 

19 



CALL ROOT (CON. W. M, 0, 
IF (KER) 22.25,19 
00 22 J-1,KER 



ABS( 1 ) - X ♦ 1.0 

0R0( 1 ) • Y ♦ (D(J )♦!. 01/2.0 

IF (NR) U5,U7,u$ 
hS 00 *6 I-1.NR 

A • ABSt 1 ) - RX(l ) 

6 • 0R0( 1 ) - RYtl ) 

IF (A«A ♦ B«B - C ) 22.22.U6 
*6 CONTINUE 
U7 CALL GRAO (PT.ABSt 1 ).g*Qtl).OX,OY, QUAO, OELX, OELY, OS, H,0) 

IF (QUAO - l.E-10) 22,22,26 

20 IF (OELX) U8.U8.21 

21 OS » -OS 

OELY* OELY - OY - OY 
OELX- OELX - OX - OX 
IF (OELX) HB.U8,22 
U8 ABS(2) • ABS( 1) •> OELX 
OR0(2) - OR0( 1) ♦ OELY 
CALL CUT (NR, A8S,0R0,RX,RY,0S,PT,X,Y,H) 



22 CONTINUE 
3 
I • 1-26 
DO 23 J-T, 12 



DO 23 1-1.19 
DELI • 1-26 



OELJ • JO 
RAD - SCRTF(OELI«OELI ♦ CELJ-DELJ) 
ANGLE • HKI I • J ) •! Q 
ANGLE • ANGLE/57. 2957 7951 
COS - CCSF(ANGLE) 
SIN • SINF(ANGLE) 

COSK(I,J) - -(OELl«COS ♦ DELJ»SIN)/RAO 
23 SINK(I.J) - (0€LX»S;n - D£LJ»COS)/RAO 
IXP • IXKIN ♦ 1 
IYP - IY«IN ♦ 1 
DO 2% IY • !YP, IYMAX 
Y - IY 

DO 2% IX ■ IXP,IX*AX 
X - I X 
COS - COSK(IX.IY) 



SIN 

ABSI 1 



SIN 



28 
2U 
29 



30 



32 
25 



7 
8 



10 



0R0( 
ABS( 
0RC( 
APS( 
CRD( 
ABM 
CR0( 
A8S( 
C«D( 

00 is 

A8S( I 
0R0( I 
CALL 
CCNTI 
READ 
?OR*A 
READ 
00 30 
ABSI 1 
ORD( I 
CALL 
REAO 
REAO 
00 32 
ABS( I 
OROd 
CALL 
STOP 
END 
SUBRO 
DIHEN 
P( 1 ) 
P(2) 
LX » 
DO 1 
P(2M 
P(6*I 
P(U) 
PI 12) 
CO : J( 1 
C0N12 
C0NI3 
CONIU 
CON! 1 
C0NI2 
CON! 
CON! 
CON! 
CON( 
CON! 
CON! 
CON! 1 
C0NI2 
C0NI3 
CONIU 
ENO 
SUPRC 
OIMEN 
FOURT 
IF (C 
IC 
(C 
(F 
1 ) 
2) 



) 
) 

) - 
) • 
) • 
) • 
) • 
) • 
) - 
) - 

I- 
) - 
) • 
CRAM 
NUE 

?'l 
"!. 

) ■ 

) • 
CRAM 

"I- 

) ■ 

) • 

CRAM 

25 

UTIN 
SIGN 
« HT 
• HT 
IX-2 
1-1 



* ( 
x 

Y ♦ 
X ♦ 

Y ♦ 
ABSI 
CRQI 
ABSI 
OROI 
ABSI 
ORDi 

0.*» 



IX. IY) 

♦ C C S/ 3 . 
SIN/3.0 
x - ABSI 
Y - ORC( 



I ) 
1 ) 



21 

:> 
i) 
i) 
:> 
:) 



M.73; 
1 1.73; 



♦ SJN/9.0 
- COS/9.0 



050807«COS - 
C50807-SIN ♦ 



SINJ/13.0 
CCS)/18.0 



ABS( I ) - *»,0 



.9»CR0( I ) - 3.2 

(S.A8S. OR 0,2,0, LA, IT, !.,!.. 0,0, 2, 2,8,8, 0,LAST) 



I AB 

U(T 

IORO 

1.52 

0.8» 

O.P» 

(52 
( ABS 
IORO 
1.52 
0.8* 
0.8« 

(52 



(I). 1-1.52) 

F6.3/) 

(I). 1-1.52) 

ABSI I ) - W.O 
0R0( I ) - 3.2 

,4BS,ORO,2.0,LA,IT,l.,l.,0,0,2,2.3,a,0,LAST) 
(I). I-T.5J) 
(I). 1-1.52) 

ABSI I) - h.Q 
ORDI I ) - 3.2 
,AaS,0R0,3,0.LA,IT,l.,l.,0,0,2,2,8,8,0,LAST) 



E COEF I 

HT( 19. 

I IX, IY-1 

( IX^l.IY 



HT. IX, 
12). PCI 
) 
-1) 



i! 



V. P. CON) 
)t CCNIU.U) 




IF 
IF 
IF 

01 
C( 



LT 
SI 
H 

CN 
CN 
CN 
CU 



INE ROOT (CCN, KS, H, D. K£R) 
ON CON(U.U), 0(3). X(3) 



CONIU.KS) - 16.0-H 
12.2.12 



(l.KS) ) 
(2.KS) ) 
(3.KS) ) 
RTH) 
-1.0 
1.0 
2 



6,3. 6 

lit*. 11 
9.5. 9 



KER 

RETl'RN 

B - 0.5«C 

C « FCURT 

OUAO - B» 

RAD ■ SC 

IF (OUAO) 

IF IRAO - 

X( 1 ) - -B 

KER « 1 

CO TO 28 

KER • 

RETURN 

X( 1 ) - -8 ♦ RAO 



)/C0N(2,KS) 
, KS) 



0N(3,K 
H/CONI 
8 - C 

RTF! ABSFIQUAO)) 
7,8.10 
1.E-5) 8.8.9 



i 



2: 

2: 

2» 

2! 

2< 



1 

2 

5 

6 

7 

8 

9 



- 16 - 



11 



12 



13 



15 



16 



17 
18 
19 

22 
2U 
25 



26 
27 



26 



29 

30 



I 



X(2> • -b - RAO 
KER • 2 
GO TO 28 

X( 1 ) ■ -F0URTH/C3NI3.KS) 
KER ■ 1 
CO TC 28 

A • CCNC,KS)/CON( ItKS) 
CCN(i,*S)/CON( 1,k£) 
FCURTh/CON( 1,KS) 
A«A/9.0 - 6/3.0 
A«e/6.0 - C/2.0 - A»A»A/27.0 
• SCRTF(ABSF( ?)) 
(P) li, 16,18 

3« A5INM(0/(RA0»RA0»RA0) 1/3.0 
) • 2.0»RA0«SINH(PMI3) - A/5,0 

(1.732050807»RAO«COSH(PH13) - l.E-5) 15 t 15,28 
) • -<X( l)*A)/2.0 

■ 2 
TO 28 
) ■ 72. 0»ABSF(Q))»«. 33333333333333 - A/3.0 

(Q* 17.28,28 

) • -xt 1 ) - A/ 1.5 

TO 28 

- Q/(RAO«RAO*RA0) 
AeSF( ARC) - 1 .0) 27.25, 19 

- ACOSH(ABSF(ARG» 1/3.0 



P 

• 
RAO 

PHI 
X( 1 
KER 
IF 
X(2 
KER 
CO 
X( 1 
KER 
IF 
X( 1 
GO 
ARG 
IF 
PHI 
IF 
X( 1 
KER 
IF 
X( 1 
GO 
X( 1 
X(2 
KER 
IF 
X(2 
GO 
PHI 
SIN 
Z 
X( 1 
X(2 
X(3 
KER 
JAR 
00 
IF 
JAR 
0( J 
CCN 
KER 
ENO 
SUP 
OIM 
XX 
Y* 
XP1 
XMl 
YP1 
YM1 
TX 
TY 
CIR 
X2M 
Y2M 
SXP 
SXH 
S>P 
SYH 
QXP 
QXM 
OYP 
QYM 
01 
02 
03 
QU 
PI 



S 

(C) 21.2U.2U 

) - -2.6«RA0»C0SH(PHI31 - A/3.0 

» 1 

( 1 .732050807»RA0»SINH(PHI3) - l.E-5) 15,15.28 
) » 2.0«RAD»COSH(PHI3) - A/3.0 
TO 22 

> • 2. 0«A0SF(3)". 33333333333333 - A/3.0 
) ■ -<X( l)*A)/2.0 

« 2 
(Q) 26.28,28 
) ■ -x(2) - A/ 1.5 
TO 17 
3 - ACOSFf ARG) /3.0 

« 1.732C50807»RA0^SINF<PHI3I 

■ RA0«C0SF(PHI3) 
) - Z ♦ Z - A/3.0 
) « -Z ♦ SIN - A/3.0 
) * X(2) - SIM - SIN 

- 3 
» 
30 IM,KER 
(ABSF(X( I) ) - 1.0) 29,29,30 

» JAR ♦ 1 
AR ) * XII) 
TINUE 
■ JAR 



RO 
EN 



UTINE 
SION 
2.0MX 
2.0MY 
XX ♦ 1 
XX - 1 
YY ♦ 1 
YY - 1 

2.0»xx 

2.0«YY 

2.0»(X 
XPl^X 
YP1 «Y 
6.C»X 
SXP - 
6.0Y 
SYP - 
18. 0« 
OXP - 
18. 0« 
CYP - 
X2H»P 
X2M»P 
X2M«P 
X2H«P 

(TX»P( 



GRAO (P, X, Y, OX, DY, QUAO, OELX, OELY, OS 
P( 12) 

- INTF(X)) - 1.0 

- INTF(Y) ) - 1.0 

.0 
.0 
.0 
.0 



X«XX-*YY«YY) - 10.0 

HI 

HI 

X 

u 

Y 

U 
XX ♦ 2 .0 

u. 

YY ♦ 2.0 
U 

( 1 

( 
( 



KC) 



♦ 2. 

.0 

♦ 2. 
.0 

♦ 2 
.0 

♦ 2 
.0 
0) ♦ 



7) 

3) < 
6) < 
0) - 



Y2H«P(12) - 
Y2M»P( 11)- 
Y2H»P( 1) - 
Y2H»P( 2) - 
SXH.P(9) )»XP1 



<CIR-TX-TY)«P(9) 
<CIR*TX-TY)«P(8) 
(CIR*TX*TY)«P(U) 
(CIR-TX^TY)«P(5) 

- 17 - 



1 



92 • (TX»P( 7) - SXP»P(8> )»X*1 

P3 • (TX«P( 31 - SXP«P(U) )«XN1 

PU a (TX»P( 6) - SXM»P(5) )»XP1 58 

P5 ■ (TY»P(12) - SVM«P(9) )»YP1 29 

P6 • (TY«P( 2) - SYP»P(5) )«YM1 

P7 • (TY»P( 1) - SYP«P(U) )»YH1 31 

Pft • (TY«P|I!) - SYM»P(8) 1«YP1 ll 

hx ■ (YP1«(Q1-Q?*P1-P2) ♦ YMJ»(Q3-CU*P3-P«i))/16.0 33 

HY « (XPW(Q1-Q**P5-P61 ♦ XH1«(Q3-Q2*P7-P81)/16.0 3% 

OUAO ■ MX«HX ♦ HY»HY 35 

IF (QUAO - 1.1-10) 3,3,2 ; 6 

2 IF (KC) 5.U. S 37 

k HXX • ( VP 1«[§XP»P ( 10)-QXP«P(9)-SXM«P( 7)«QKP«P(f ) ) ♦ 38 

♦ YMU(SXM»P( 3)-QXM«P(U)-$XP»PC6>*QXP«P(5) 11/8.0 39 
MYY ■ (XP1»(SYP»P( 12)-QYP«P( *)-SYM*P( 2)*QYH«P(5) 1 ♦ *»0 

♦ XMU(SYM»P( 1)-QYN»P(H)-$YP»P( lb*QYP«R(8)) 1/8.0 *T 

HXY - (Pl-PU*P3-P24PS-P8*P7-P6*Ql-Q2*Q3-QM/%,0 ** 

RAO - SQRTF(QUAO) %3 

DX « DS»HY/*AD %% 

OY ■ -OS«HX/RAO *5 

SEC « 0.5*(HXX»DX»DX ♦ HXY«OX«OY ♦ HYY»OY»OY) /QUAO *6 

D6LX« DX - HX»SEC %7 

DELY- OY - HY«SEC **8 

RETURN *9 
5 OH - H - ( YP1«(XP1«Q1-XH1«Q21 ♦ YMU(X«l»Q3-XP1»0%!)/32.0 
OX ■ CH»HX/QUAO 
OY ■ CH«HY/QUAO 
OELX ■ OH«OH/QUAO 

3 ENO 
SUBROUTINE OUT <NR . ABS.ORO.RX .RY, OS, PT. X.Y.HI 
01 HENS I ON ABS ( 900 ), 0R0( 9001 ,RX(%81,RYUdl 9 PT( 121 

C « 2.0«DS«DS 1.5 
NPT - 1 
NR « NR ♦ 1 

RX(NR) • ABS( 1) % 

RY(NR) - OR0( 1) S 
GO TO 1 

U IF (UC-NPT) 20,20, 2U 

2U CALL GRAO (PT,A8S(NPT),QR0(NPT1,DX,0V, OUAO, OELX, OELY, OS, H, 01 8 

IF (OUAO - 1.E-10 ) 19,19.5 9 

5 AB$(NPT*1) • ABS(NPT) ♦ OELX U 

0RD(NPT4l) ■ ORO(NPT) ♦ OELY 15 
CALL GRA0(PT,ABS(NPT*1) t 0R0<NPTO 1 ,OX,OY, OUAO, OELX, OELY,0$,H, 1 1 



IF (QUAO - 1.E-10 J 19.19,6 

15,25,19 



i 



* 



6 IF (OELX - C) 25,25,19 18 

25 NPT » NPT ♦ 1 19 
ABS(NPT) « ABS(NPT) ♦ OX 20 
ORO(NPT) * ORD(NPT) ♦ OY 20.1 
NPT1« NPT % 21 
IF (ORD(NPT)-Y) 7,10,8 22 

7 NPT1« NPT1 - 1 23 
GO TO 10 2% 

8 IF (Y*1.0-OR0(NPT ) 1 9,10,10 25 

9 NPT1- NPT1 - 1 

10 IF ( AeS(NPTIl-X) 11,U,12 

11 NPT1- NPT1 - 1 |8 
GO TO 15 29 

12 IF (X*1.0-ABS(NPT1 )) 11,U,U 
1U IF (NPT-NPT1) 15,U,15 
15 NPT « NPT1 

19 IF (NPT-1) 26,26,20 

20 NR ■ NR ♦ 1 3* 
RX(NR) ■ ABS(NPT1 35 
RY(NR) - QRD(NPT) 36 
00 21 1*1, NPT 37 
ABS(I) - 0.8»ABS(I) - U.O J8 

21 ORO(I) - 0.8«0R0(I) - 3.2 39 
LA » UH 

CALL CRAW (NPT,A8S,0R0,2,0,LA,7,1.,1.,0,0,2,2,8,8,0,LAST) 
RETURN 

26 NR - NR - 1 
ENO 
ENO 

10 21 3 6 13 h 10 {EROTH 



- 18 - 



J! 












0F1 

END 



1 



10 !» m 

I/ 8 



11 



1-25 !| 1-4 iv e 






I 

1-1 t «ii "■ 
'l^S 11-4 »V6 



UNCLASSIF'^ 



DUDLEY KNOX LIBRARY - RESEARCH REPORTS 



5 6853 01068951 6 










DEFENSE TECHNICAL INFORMATION CENTER 
CAMERON STATION 
ALEXANDRIA, VIRGINIA 22304-6145 



OFFICIAL BUSINESS - PENALTY FOR PRIVATE USE $300 
POSTMASTER: DO NOT FORWARD 



Au Number 
604353 



Class 
U 



Pages 




iype Copy 
H 



Received Date: 14-SEP-9S 

Requested By : JODY OTTILIGE 
ATTENTION : DR M RENEKER 

NAVAL POSTGRADUATE SCH 
MONTEREY, CA 93943-5000 




CODE 524 -LIB DIR 



UNCLASSIFIED